Draft version October 26, 2012 

Preprint typeset using L£-T^X style cmulatcapj v. 03/07/07 



LOCAL STUDY OF ACCRETION DISKS WITH A STRONG VERTICAL MAGNETIC FIELD: 
MAGNETOROTATIONAL INSTABILITY AND DISK OUTFLOW 

Xue-Ning Bai 1 ' 2 ' 3 & James M. Stone 1 

Draft version October 26, 2012 

ABSTRACT 

We perform three-dimensional vertically-stratified local shearing-box ideal MHD simulations of the 
magnctorotational instability (MRI) that include a net vertical magnetic flux, which is characterized 
by midplane plasma (3 (ratio of gas to magnetic pressure). We have considered fio = 10 2 , 10 3 and 10 4 
and in the first two cases the most unstable linear MRI modes are well resolved in the simulations. We 
find that the behavior of the MRI turbulence strongly depends on [3q : The radial transport of angular 
momentum increases with net vertical flux, achieving a ~ 0.08 for j3 = 10 4 and a > 1.0 for fto = 100, 
where a is the height-integrated and mass- weighted Shakura-Sunyaev parameter. A critical value lies 
at /3o ~ 10 3 : For /?o > 10 3 , the disk consists of a gas pressure dominated midplane and a magnetically 
dominated corona. The turbulent strength increases with net flux, and angular momentum transport 
is dominated by turbulent fluctuations. The magnetic dynamo that leads to cyclic flips of large- 
scale fields still exists, but becomes more sporadic as net flux increases. For /3q < 10 3 , the entire 
disk becomes magnetic dominated. The turbulent strength saturates, and the magnetic dynamo is 
fully quenched. Stronger large-scale fields are generated with increasing net flux, which dominates 
angular momentum transport. A strong outflow is launched from the disk by the magnetocentrifugal 
mechanism, and the mass flux increases linearly with net vertical flux and shows sign of saturation at 
Po i$ 10 2 . However, the outflow is unlikely to be directly connected to a global wind: for (3q > 10 3 , 
the large-scale field has no permanent bending direction due to dynamo activities, while for /3q < 10 3 , 
the outflows from the top and bottom sides of the disk bend towards opposite directions, inconsistent 
with a physical disk wind geometry. Global simulations are needed to address the fate of the outflow. 

Subject headings: accretion, accretion disks — instabilities — magnetohydrodynamics — methods: 
numerical — turbulence 



1. INTRODUCTION 

The magne t orotat ional instability (MRI, 

iBalbus fc Hawlevl 1 1 9 9 If) is believed to be the pri- 
mary mechanism for driving accretion in a wide range 
of astrophysical disks by generating turbulence to 
provide outward angular momentum transport. The 
robustness of the MRI and its physical properties have 
been widely stud ied by means of both local shearing box 
simul ations (e.g.. lHawlev et af1ll995t [Brandenburg et all 
19951) and global simu lations (e.g^ lHawlevI 120011 : 
Fromang fc Nelson! 120061 ) . Local shearing-box simula- 
tions have the advantages of being able to achieve high 
resolution at modest computational cost, and are the 
primary way for exploring the detailed physics of the 
MRI. 

Four types of shearing-box MRI simulations exist, 
depending on whether vertical gravity from the cen- 
tral object (which gives density stratification) is in- 
cluded and whether the simulations include net ver- 
tical magnetic flux. While unstratified MRI simula- 
tions provide important benchmarks on the physical 
prope rties of the MRI t urbulence (e.g., lHawley et all 
119961: ISano et all 12004 iFromang fe Papaloizoul l2007t 
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IFromang et all l2007t iLesur fe Longarettil |2007( ). real 
disks involve vertical density stratification. So far most 
vertically stratified shearing-box simulations focus on a 
magnetic field c onfiguration with zero net vertical mag- 
netic flux (e.g.. IStone et all 119961: IMiller fe Stone! 120001: 
Zieeler fe Riidigerl 120001 : iHirose et al l 120061 : iDavis et al 
20101: iShi et all 120101: iFlaig et all 12011 iSimon et al 



20121 ). However, the inclusion of net vertical magnetic 



field in stratified shearing-box simulations, which is prob- 
ably closer to reality, has only been ra rely explored , most 
likely due to numer ical difficulties (|Stone et all 119961 : 
IMiller fc Stond [20001 ): either the numerical code crashes 
or the disk is violently disrupted in a few orbits. 

More and more attention has been drawn to this least 
explored regime of stratified shearing-box with net ver- 
tical flux in the recent years. This is partly due to the 
improvement of the numerical schemes, but more impor- 
tantly, the physical significance of the MRI in this regime 
and its intimate connections to astrophysical applications 
make it deserve more detailed study. In particular, the 
presence of jets and outflows in a wide range of accretion 
systems such as protostars fe.g.. iReipurth fc Ballvll2001l: 
Cabritll2007l) X-ray binaries (e.g.. iMirabel fc Rodriguez! 



19991: iFenderl l2006f) and active galactic nuclei (e.g. 
Begelman et allll984l : lHarris fc Krawczvn ski 2006) is in- 
dicative of the presence of large-scale (poloidal) mag- 
netic field threading through the disk, which is es- 
sential for jet/outflow acceleration and collimation via 
the magnetocentrifugal mechanism (|Blandford fc Pavnd 
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fl98l iPudritz fe NormanlH983l) 

Rec ently, ISuzuki fe Inutsukal ()2009l ) and ISuzuki et al.l 
()2010() successfully conducted stratified shearing-box 
simulations that include net vertical magnetic held with 
an outflow boundary condition using characteristic de- 
composition. Their net vertical field, characterized by /3o 
(ratio of gas to magnetic pressure at midplanc), is rather 
weak (/3q > 10 4 ). They reported the launching of disk 
outflow and found that the averaged mass outflow rate 
increases linearly with the magnetic pressure of the net 
vertical field. It was speculated that the outflow serves as 
the base of a Blandford-Payne type magneto-centrifugal 
wind. 

In this paper, we conduct shearing-box simulations of 
the MRI that include a relatively strong net vertical mag- 
netic flux, with 10 4 > f3 > 100 which has largely been 
unexplored in the literature. The main motivation of our 
work is two-fold as we elaborate below. 

First, we aim at studying the properties of the MRI 
turbulence in more realistic magnetic field geometry 
and their dependence on the net poloidal magnetic flux 
threading through accretion disks. In particular, how 
the Shakura-Sunyaev a parameter, which is of the most 
common interest, depends on the net magnetic flux. Due 
to the vast dominance of zero net vertical flux strati- 
fied shearing-box simulations, the a parameter has been 
taken for granted to be of the order 0.01, while the a 
value inferred from observations of fully ionized disks is 
at least an order of magn itude higher (e.g., see discus- 
sions by iKing et ail 120071 and references therein) . Un- 
strat i fed shearing-box s imulations (e.g., lHawlev et al.l 
Il995t iBai fc Stond 1201 ll ) unambiguously found that as 
long as the most unstable linear MRI mode fits into the 
simulation domain, the turbulent stress increases roughly 
linearly with 1//3q. It is natural to expect the same trend 
in stratified shearing-box simulations, hence net verti- 
cal magnetic flux is likely to be a crucial ingredient in 
real accretion disks. Moreover, as in the zero net-flux 
cases, we expect buoyancy and open vertical boundaries 
in our stratified simulations to give rise to novel features 
on the properties of the MRI turbulence, particularly 
the MRI dynamo and the generation o f large-scale fields 
dLesur fc Qgilvid IxM iGressell [20101 : iGuan fc Gainlmel 
[2Q11D. 

Second, our study will address the potential connec- 
tion bet ween the MRI and the ma gneto-centrifugal wind 
(MCW, iBlandford fe Pavnd Il98l . The MCW extracts 
both mass and angular momentum from accretion disks, 
and is a competing mechanism for driving disk accretion 
and evolution. While the wind acceleration and colli- 
mation of the MCW is largely known, the wind launch- 
ing process which loads mass from the disk onto the 
wind is still not well understood. Most semi-analytical 
works and numerical simulations either assume a razor- 
thin disk treated as a bound ary condition with arti- 
ficial mass injection (e.g . . IlI [19951 : lOuved fc Pudritzl 
Il997t iKrasnopolskv et all I1999I ). or proceed with un- 
resolved disk by adopting some forms of artificial dif- 
fusion, (e.g iKato et al.l 120021: ICasse fc Keppensl 120021: 
IZanni et alj |2007D . Such diffusion, which presumably 
arises from the MRI, is necessary to allow the gas to slide 
through the disk that avoids rapid accumulation of mag- 



netic flux to maintain steady state accretion. However, 
the microphysics (e.g., the MRI) within the accretion 
disk, which is crucial for wind launching, is not mod- 
eled directly. In fact, it is numerically very challenging 
to properly resolve the MRI turbulence in global wind 
simulations, which requires extended three-dimensional 
computational domain with large numerical resolution. 
On the other hand, local shearing-box simulations pro- 
vide superb resolution with realistic computational cost, 
hence are a powerful tool for studying the wind launching 
process fro m first principle . 

Recently, iFromang et al.l ()2012l ) conducted a series of 
stratified shearing-box MRI simulations with /3o = 10 4 . 
After carefully examining various numerical issues, they 
reported simultaneous MRI turbulence and launching 
of an MCW, whereas angular momentum transport is 
still do minated by the MR I rather than the MCW. In 
parallel, iLesur et ail (|2012l ) conducted a series of strati- 
fied shearing-box simulations in both one-dimension and 
three-dimensions with /3q ~ 10, where the most unstable 
linear MRI modes just marginally fit into the disk. They 
found that such MRI modes directly produce magneti- 
cally driven outflows, which are time varyin g due to sec- 
ondary instabilities. Meanwhile, iMollI (|2012[ ) studied the 
wind launching process in two-dimensional shearing-box 
simulations with /3q ~ 1 — 10 (where MRI is suppressed) 
and confirmed the mass-flux instability in t he strong field 
regime found in earlier a nalytical works (|Lubow et al.l 
119941 : ICao fc SpruTtll2002f ). Our simulations fill the gap 
in parameter space explored by these authors and reveal 
the interesting transition from the MRI-dominated trans- 
port to potentially wind-dominated transport. Together 
wi th the e a rlier s emi-analytical study of wind launching 
bv I Qgilvid ()2012l ). which appl ies to the regime of / 3n 5, 1 , 
as well as the simulations by ISuzuki fc Inutsukal (|2009T) , 
which applies to /3 > 10 4 , the series of studies cover the 
complete parameter range relevant to MRI turbulence 
and wind launching in realistic accretion disks under the 
local shearing-box framework. 

This paper is organized as follows. In Section [2] we 
describe the methodology and setup of our simulations. 
We present our simulation results in Section [3l focusing 
on the properties of the MRI and angular momentum 
transport, and Section IH focusing on the dynamics of 
the outflow and its connection to global disk winds. Each 
result section is broken into three subsections, devoting 
to one aspect of the results. Implications of our results 
for global disk evolution arc discussed in Scction[5]bcforc 
we conclude. 

2. SIMULATION SETUP 

We use the Athena magnetohydrodynamic (MHD) 
code (jStone et al.ll2008D and perform three-dimensional 
(3D) local MHD simul ations of gas dynamics using th e 
shearing-box approach (jGoldreich fc Lvnden-Belllll965| ). 
Ideal MHD equations are written in a Cartesian coordi- 
nate system in the corotating frame at a fiducial radius 
R with Keplerian frequency £1, with e x ,e y ,e z denoting 
unit vectors pointing to the radial, azimuthal and vertical 
directions respectively, and ft is along the e z direction. 
The MHD equations read 

g+V-(/m) = 0, (1) 
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dpu 



V • (pu T u + T) = p 



2uxfl + 3fl 2 xe x - Vl 2 ze z 



dB 

— =Vx(uxB), 

where T is the total stress tensor 

T= (P + B 2 /2) I — B T B 



(2) 
(3) 

(4) 



TABLE 1 

Summary of Simulation Runs. 



p is the gas density, u is the gas velocity. We use an 
isothermal equation of state P = pc 2 . Note that the 
unit for magnetic field is such that the magnetic pres- 



sure equals B /2, which avoids the extra factor of V47r. 
The vertical gravity from the central object (— Q 2 z) is in- 
cluded in the momentum equation hence our simulations 
have vertical density stratification. 

The shearing-box source terms (Coriolis force and 
tidal gravity) have been readily implemented in Athena 
(jStone fe Gardine"rll2010[ ). An orbital advection scheme 
is adopted where the system is split into an advec- 
tive part that corresponds to the background shear flow 
— 3flx/2e y , and the other part that only evolve the ve- 
locity fluctuations v: 



v = u 



2 v 



(5) 



Most of our analysis deals with v, except in the study of 
conservation laws along field lines, where u will also be 
used (see Section |4.1.2[) . 

We use standard s hearing-box bounda ry conditions in 
the radial direction (jHawlev et alJll9 95). In the verti- 
cal direction, we adopt a simple zero-gradient outflow 
boundary condition for velocities and magnetic fields, 
with density extrapolated assum ing vertical hydrostatic 
equilibrium (jSimon et al.l l20lil ). In addition, vertical 
velocity in the ghost cells is set to zero in the case 
of inflow. It was pointed out that strict zero-gradient 
boundary conditio n would prevent MHD-driven outflow 
(|Lesur et al.ll2012"l ) , while we find that the outflow can be 
launched once the density extrapolation is implemented 
in the ghost zone. 

It is well known that in the presence of pure verti- 
cal background magnetic field, the linear MRI modes 
consist of counter-motions in the radial direction that 
alternate along the vertical direction (i.e., the channel 
flow) . The channel flow is an exact solu tion of the MHD 
equat ions even in the non-linear regime (|Goodman fe Xul 
1994). It can achieve very large amplitude before 
broken up by parasitic in stabilities (lLatter et al.l l2009t 
iPessah fe Goodmanll2009t ). Recently. lLatter et al.l(|2010H 
studied the MRI modes in the presence of vertical strati- 
fication and their parasitic instabilities. They found that 
the channel modes strongly amplify the magnetic fields 
to thermal strength before they can be destroyed. On 
the simulation side, the channel flows in the initial stage 
of the simulations not only place enormous demands on 
numerical codes, but also give rise to er uptive behavior 
that depletes the disk in just a few orbits ([Miller fc Stond 
lloooh . in view of the transient nature of the channel 
flow, as well as its large numerical demand, we initialize 
our simulations in a way that strong channel flow can be 
avoided, as described below. 



Run 


Po 


Domain Size 


Resolution 


Duration 1 


B2 


10 2 


4H x 8H x 12H 


96 x 96 x 288 


960C- 1 


B3 


10 3 


4H x SH x 12H 


96 x 96 x 288 


960Q" 1 


B4 


10 4 


4H x SH x 12H 


96 x 96 x 288 


960Q" 1 


B3-hr 


10 3 


3H x GH x 12H 


144 x 144 x 576 


720ft- 1 



■ duration of the run with net flux fully added. 

Fiducially, we use a box size of (4i/, 8H, 12H) in 
(x, y, z) resolved by (96, 96, 288) grid points, where H = 
c s /n is the thermal scale height. The relatively large 
horizontal box size is necessary to properly capture the 
mesoscale stru ctures of the MRI turbulence such a s 
the zonal flow ([Johansen et al.1120091 ; iSimon et al.ll2012l) . 
The initial density profile is taken to be Gaussian p = 
Po exp (— z 2 /2H 2 ). We use natural unit in our simula- 
tions, with c s = 1, CI = 1, H = 1 and po = 1. We 
initialize our simulations with zero net vertical magnetic 
flux configuration, i.e., a weak vertical magnetic field that 
varies sinusoidally in the radial (x) direction. We run the 
zero net-flux simulation to t = 120f2 -1 when turbulence 
is fully developed, then we start to add net vertical mag- 
netic flux, which is applied uniformly in the grid (hence it 
does not introduce any divergence error) and the amount 
added is proportional to the time step dt. This process 
continues to t = 240f2 _1 when the net vertical magnetic 
field strength reaches the desired level Bq, characterized 
by midplanc plasma /? = 2p c 2 /B 2 . The simulations are 
run for about another 150 orbits to t = 1200r2 — 1 . 

For our simulations, we consider /3q = 10 2 , 10 3 and 10 4 , 
corresponding to runs B2, B3 and B4. Ignoring vertical 
stratification, the most most unstable MRI wavelength 
at di sk midplane is about 9.18(5q 1 ^ 2 H (jHawlev et al.1 
19951). Including v ertical stratification, calculations by 
Latter et alJ (|2010l ) suggest resolution of about 25 cells 
per H for /3 = 100, and 50 cells per H for /3 = 10 3 , 
to properly resolve the MRI channel modes. Our fidu- 
cial resolution meets this criterion for /3q = 100. For 
/3q = 10 3 , we conduct an additional simulation (run B3- 
hr) with twice the resolution which then does properly 
resolve the most unstable linear modes. For this run, we 
initialize the simulation with full a net vertical magnetic 
flux, together with a sinusoidally varying vertical field 
component so as to avoid an excessively strong initial 
channel flow. The list of simulation runs is provided in 
Table [TJ 

A density floor of p F ioor = 10~ 4 ,3 x 10~ 5 and 10~ 5 
are applied for simulations with /3q = 10 2 ,10 3 and 10 4 
respectively to avoid numerical difficulties at magneti- 
cally dominated (low plasma /3) regions. We will see in 
Figure [3] that horizontally averaged densities in the satu- 
rated states of all our simulations are at least two orders 
of magnitude larger than the density floor. In addition, 
we apply another density floor that is ten times smaller 
than pfi oor to all the left and right states inside the MHD 
integrator (we use the CT U integrator with third order 
spatial reconstruction, see lStone et al.ll2008l for details). 
This procedure is essential and greatly improves the ro- 
bustness of the Athena MHD code. 
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Fig. 1. — One snapshot from our simulation run B3. The surface 
plot shows the isosurface of gas density, increasing logarithmically 
from blue with p = 0.1 (cutoff) to red with p = 2.0. The stream- 
lines indicate gas velocity, increasing linearly from white with v = 
to dark red with v = 5c s . 

Our simulations with /3q — 10 2 and 10 3 occasion- 
ally produce extremely strong magnetic field in very lo- 
calized regions (a few grid cells in the entire compu- 
tational domain), which severely limits the simulation 
time step. In order to proceed, we have included some 
hyper-resistivity r\ oc J 2 in regions with excessively large 
J, where J is the current density. More specifically, 
we define jd = |V x B\/y/p, and set Ohmic resistivity 
f] = 0.1(c^/fi) • (j d /10 3 fi) 2 when j d > 10 2 O and 17 = 
otherwise. We have examined that for the f3 = 10 2 and 
10 3 cases, hyper-resistivity is applied to at most 1% of 
all grid cells (most of the time less than 0.1%) and does 
not affect the general properties of the MM turbulence. 

Our simulations always launch a strong outflow (see 
Figure |9]) that would drain the gas in the simulation box 
in tens to hundreds of orbits. In order for the system 
reach a steady state which allows us to characterize the 
dynamical structure of the turbulent disk, we enforce the 
total mass contained in the simulation domain to be con- 
stant by multiplying the density of all grid cells by a 
common factor after each time step. The momentum re- 
mains fixed in each grid cell. This approach is similar to 
lOgilviel (|2012f) who considered much stronger magnetic 
fields that would launch even stronger disk outflow. As 
will be discussed briefly in Section 15.2.31 such mass ad- 
dition does not significantly affect the energetics of the 
system for all our simulation runs. 

In Figure [1] we show a snapshot of our simulation run 
B3 in the saturated state of the MRI turbulence. We 
see clearly the large density fluctuations as a result of 
the vigorous MRI turbulence, and a substantial fraction 
of the gas is in supersonic motion. Meanwhile, there is 



a systematic vertical velocity at the vertical boundaries, 
indicating a strong outflow. In the following two sections, 
we discuss in detail the properties of the MRI turbulence 
and the outflow respectively. 

3. SIMULATION RESULTS: MRI TURBULENCE 

One of the main diagnostics for the MRI turbulence is 
the rcj) component of the stress tensor, consisting of the 
sum of Reynolds and Maxwell stresses, 



T 



Roy 



pv x v y 



T, 



Max 



-B T B n 



(6) 



and they characterize the rate of radial angular momen- 
tum transport. We are interested in the horizontally av- 
eraged vertical profiles (notated with an over bar, here 
T r( p) of the stresses, and they are usually normalized 
by poc 2 .. Another useful measure is the mass weighted 
average of the stresses, no rmalized to the gas pressure 
(jShakura fc Sunvaevl 11971 



/ T r< t,(z)dz 
J p{z)c 2 dz 



(7) 



where a denotes either Reynolds or Maxwell stresses 
(with subscript), or their sum (without subscript), which 
characterizes the total rate of radial angular momentum 
transport in the disk. 

We show the evolution of aR oy and ajyiax in Figure 
[2j The total stress in the initial zero net-flux simula- 
tion saturates at (aR ey ) ps 5.0 x 10~ 3 and (aMax) ~ 
2.0 x 10~ 2 respectively, consistent with previous works 
(e.g., ISimon et alj|2012f ). The stresses increase steadily 
as we gradually add net vertical magnetic flux from 
t = 120Q -1 , and saturate as full net flux is in place 
at t = 240f2~ 1 . We wait for about another 20 orbits 
and take the data from t = 360f2 _1 to the end of the 
simulations at t = 1200S1 - 1 for all the time averages in 
this section, which amounts to about 130 orbits. For 
our high-resolution run B3-hr, we consider data from 
t = 180fi — 1 to the end of the simulation at t = 720f2 _1 , 
which amounts to about 85 orbits. Throughout this pa- 
per, we use angle bracket (■) to denote time average. 

In the saturated state, a common indicator for proper 
resolution of the MRI turbulence and convergence is the 
quality factor, defined as 



where 



\vAi\ = yJB]/p, (8) 



which is the ratio of a characteristic MRI wavelength to 
the gr id scale. He r e j de notes either the y or the z dimen- 
sion. ISano et ail (|2004f) found that Q z > 6 is required 
for the MRI turbulence to be well resolved. lHawlev et al.l 
(|2011f) further suggested that Q z > 10 and Q y > 20 suf- 
fice for properly resolving the MRI turbulence. Although 
these results arc mainly based on cither unstratificd sim- 
ulations or stratified simulations with zero net vertical 
magnetic flux, we take them as useful reference. In our 
simulations, the quality factors minimize at the disk mid- 
plane, where we find (Q z ) to be about 87, 50, 12 for runs 
B2, B3 and B4 respectively, and (Q y ) is about 185, 110, 
25 for the three runs. For the high-resolution run B3-hr, 
we find (Q z ) ~ 110 and (Q y ) ~ 285. In all cases, the 
quality factors are well above the suggested value. 
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t(£T 1 ) 

Fig. 2. — Time history of mass weighted Maxwell (solid) and 
Reynolds (dashed) stresses from simulations B2, B3 and B4. They 
are restarted from the end of the zero net-flux simulation at t = 
120fi _1 . Net vertical magnetic field is added gradually to its full 
strength at t = 240H -1 . 
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3.1. Disk Structure 
3.1.1. Vertical Profiles 

In the presence of a strong net vertical magnetic flux, 
the systems all saturate into very powerful MRI turbu- 
lence that substantially changes the structure of the disk. 
In Figure [3] we show the vertical profiles of various hor- 
izontally averaged quantities. The upper two panels are 
for gas density (p) and plasma (/?} (ratio of gas to mag- 
netic pressure). We see that the density profile deviates 
substantially from Gaussian due to strong magnetic sup- 
port. Only for /3 = 10 4 , the Gaussian kernel is still 
present; while for f3 < 10 3 , the entire disk becomes mag- 
netically dominated, with plasma /3 < 1 everywhere. In 
all cases, the averaged plasma /3 saturates at about 0.1 
towards the disk surface, while in the extreme situation 
of /?o = 10 2 , plasma /3 is about 0.1 everywhere as MRI 
saturates, and the disk is dramatically puffed up by the 
magnetic pressure gradient 4 . 

In the lower two panels of Figure [3] we show the pro- 
files of horizontally averaged turbulent magnetic energy 
(-EB.turb) and kinetic energy (£k,turb)- Here we have sub- 
tracted contributions from horizontally averaged mean 
fields and mean velocities at each vertical layer in the 
calculations: i?K,turb = E K — p(v) 2 /2, Es.turb = E B — 
(B) 2 jl. From the figure we see that for both magnetic 
and kinetic components, turbulent energy increases sub- 
stantially by about one order of magnitude as (3o changes 
from 10 4 to 10 3 . As net vertical field increases further, to 
/3q = 10 2 , however, turbulent energy roughly stays at the 
same level as the /3 = 10 3 case. This indicates that the strength of the MRI turbulence saturates as (3 < 10 3 . 

On the other hand, the amplitude of the mean flow and 
mean magnetic field continue to increase with net ver- 
tical magnetic field, and their contributions to the total 
energy gradually dominate that from turbulent energy. 
We will discuss this fact further in Section [XU 

For /?o < 10 3 , turbulence in all regions of the disk 
become trans-sonic or super-sonic, as one compares the 
profiles (£k,turb) and (p). This result implies strong com- 
pression and rarefaction, hence large density variations. 
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Fig. 3. — Vertical profiles of (from top to bottom) gas density (p), 
plasma (/3), turbulent magnetic energy (Bg turb) an< l turbulent ki- 
netic energy (Sk turb)- I n each panel, we show the profiles for runs 



B2 (ft) = 10 2 , solid), B3, 



lCP dashed) and B4 



10 4 



dash-dotted) respectively, and the thin dashed lines correspond to 
run B3-hr (/3q = 10 3 at high resolution). 



4 We note that ILesur et all (20121 ) and Moll ( 20H) observed 
strong compression of the disk in their simulations where stronger 
vertical field (/9o < 10) is used and reflection symmetry about the 
midplane is assumed. The compression is due to the large poloidal 
field curvature around the disk midplane. We do not find any disk 
compression in the case of /3q > 10 2 mainly because no reflection 
symmetry is enforced to build up the field curvature. In fact, our 
simulation results show opposite symmetry across the midplane 
and the poloidal field curvature near midplan e is a bout zero. See 
further discussions about symmetry in Section 14.31 
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1.2 




Fig. 4. — Vertical profiles of density fluctuations (5p 2 ) 1 / 2 (z) 
normalized to horizontally averaged density (p) for all simulation 
runs. 

Looking at the density field, we find that the density 
fluctuations in our simulations generally show elongated 
structure that is slightly tilted anti-clockwise from the 
azimuthal direction, indicatin g spiral density waves ex- 
cited by the MRI turbulence (jHeinemann fc Papaloizoul 
|2009[ ). although these structures are short lived due to 
the strong turbulence. To address such density vari- 
ations, we measure the standard deviations of density 
fluctuations bp in each horizontal layer, and normalize 
them to the horizontally averaged densities. The ob- 
tained profile of {bp 2 ) 1 / 2 / p for all simulations are shown 
in Figure |4l We sec that density variations are significant 
in all cases. For runs B3 and B4, the midplane density 
fluctuation increases with net magnetic flux, and reaches 
about 0.5 for /?q < 10 3 . The fluctuation level increases 
from the disk midplane, and reaches order unity at disk 
surface as the disk become more and more magnetically 
dominated. The case with (3q = 10 2 show distinctive 
features from higher (3o cases: except for a small bump 
around the disk midplane, the density fluctuation level 
stays roughly constant at about 0.5. We note that the 
turbulent kinetic energy profile in run B2 is similar to 
that in run B3, while the density near disk surface in run 
B2 is a factor of several higher than run B3, indicating 
smaller turbulent velocity, hence the profile of (bp 2 ) 1 / 2 / p 
for run B2 drops below that for run B3, which is essen- 
tially a reflection of weaker turbulence for run B2. 

3.1.2. Autocorrelation Function 

Another useful diagnostic of the MRI turbulence is 
the autocorrelation function (ACF) of kinet ic and mag- 
netic energy fluctuations , which is defined as ([Guan et al.l 
120091 ISimon"et~a^[20T2T) 

ACF(E K . B ,Ax) ^ f6P(t,x)dx3 A ' 

(9) 

where / denotes B for magnetic energy, and yfpv for 
kinetic energy, the angle bracket denotes time average. 
Here we consider a two-dimensional ACF in the horizon- 
tal plane at some fixed vertical height z, and we subtract 
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Fig. 5. — Two-dimensional autocorrelation function (ACF) in the 
horizontal plane for magnetic energy fluctuations for our simula- 
tions B4 (/3 = 10 4 ), B3, B3-hr (/3 = 10 3 ) and B2 (/3 = 10 2 ), from 
top to bottom. Left panels correspond to the ACFs for \z\ < 2.5H, 
while right panels are for \z\ > 2.5H. Note that the domain size 
for run B3-hr is slightly smaller. 
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the horizontally avcra_ged quantities in the evaluation of 
5f : Sf = f(x,y,z) — f(z). This is necessary as the mean 
flow and mean field can become particularly strong in the 
presence of net vertical magnetic field. 

In Figure El we show the ACFs for magnetic energy 
fluctuations in all our simulations B2, B3, B3-hr and B4, 
for both near the midplane region (left, for z < 2.5H) 
and the disk surface (right, 2.5H <z< 5H). The ACFs 
for the kinetic energy exhibits similar patterns, hence 
we do not involve a separate figure to show them. The 
ACF for run B4 at midplane where the net vertical mag- 
netic flux is relatively small (/3o = 10 4 ) is very similar 
to that found in zero n et-flux simulations (|Guan et al.l 
l2009t ISimon et al.l l2012h . containing a narrow and elon- 
gated centroid that is tilted with respect to the azimuthal 
axis by about (b ~ 15°. Such a tilt angle is related to the 
empirical correlation bet ween the plasma (3 and the stress 
param eter a « 1/2/3 (jGuan et alj l2009t iBai k Stone] 
l2011f ). which also applies for the midplane region of our 
run B4 as we can compare Figures [3] and [51 Moving up 
to the disk surface, we find that the centroid becomes 
broader in the upper layer, and the tilt angle also be- 
comes smaller. These features indicate longer correlation 
length in the radial direction, as well as more isotropy in 
the magnetically dominated disk corona. The ACFs in 
run B3 at the midplane has similar tilt angle as in run 
B4, typical for MRI turbulence, while at both midplane 
and corona, the peaks in the ACFs are as broad as the 
corona region of run B4, which is in line with the fact 
that the entire disk has become magnetically dominated. 

Our run B2 with /3q = 100 exhibits an extreme situa- 
tion: fluctuations are highly elongated in the azimuthal 
direction, with the measured correlation length compa- 
rable to the azimuthal size of our simulation box, which 
might suggest that our azimuthal box size is not sufficient 
to properly resolve the MRI turbulence. Nevertheless, 
the shape of the ACF indicates that the fluctuations are 
quasi-axisymmetric (which is also evident as one views 
the raw simulation data), and the radial fluctuations are 
well fitted into our simulation box. 



is also reflected in Table [2] where the Shakura-Sunyaev a 
from run B3-hr is about 25% higher. Similarly, in Figure 
[31 we see that run B3-hr leads to slightly larger density 
fluctuations. Later in Figure [9l slightly higher mass out- 
flow rate is achieved in run B3-hr. 

The systematically higher saturation amplitude of the 
MRI in the high-resolution run discussed above might 
suggest that our simulations for (3q = 10 3 have not fully 
reached numerical convergence. However, such higher 
saturation amplitude may not be due to under-resolved 
MRI turbulence, and several other factors are likely to 
play a role. While slightly smaller horizontal box size 
in run B3-hr can be a potential cause, the fact that 
the peaks of the ACF are well contained in t he box 
makes it somewhat unlikely ([Simon et al.1 l2012f ) . An- 
other possibility is that in stratified MRI simulations, 
the property of the MRI turbulence is largely affected by 
the dynamo process (see detailed study in Section [3~3)l . 
which generates strong mean toroidal field that changes 
sign over time. It appears that the instantaneous large- 
scale toroidal field generated in run B3-hr is on average 
stronger than that in run B3, as can be judged from Fig- 
ure [7j and both are above equipartition field strength. 
Such larger mean toroidal field leads to stronger mag- 
netic support to the disk, resulting in slightly different 
vertical disk structures between the two runs B3 and B3- 
hr. In this respect, one should be careful when discussing 
the issue of convergence, since comparing turbulent prop- 
erties between simulations with different disk structures 
can be misleading. Since the mechanism of the dynamo is 
far from being understood, and the dynamo in the pres- 
ence of strong net vertical flux behaves very differently 
from the case with zero net vertical flux, we defer this 
convergence issue for future study. 

In sum, although quantitative differences exist between 
runs B3 and B3-hr, the physical properties of the MRI 
turbulent disk from the two runs all agree with each 
other. We therefore consider both runs as valid repre- 
sentations of the physical system, with run B3-hr as a 
more accurate model. 



3.1.3. Numerical Convergence 

Due to the large computational cost, we have not car- 
ried out a full resolution study in this paper. How- 
ever, we do expect our simulation results to converge 
in runs B2 and B3-hr where we have adequate resolution 
to properly resolve th e linear MRI modes as suggested 
bv lLatter et al.l (l2010f). B ased on the quality factor crite- 
rion ( Hawlev et al.ll201~ll) . we expect all our simulations 
to well resolve the MRI turbulence. Below we compare 
the results between our simulation runs B3 and B3-hr 
and briefly discuss about numerical convergence. 

We see from Figure [5] that the shape of the ACFs in 
our runs B3-hr and B3 are similar, but the ACF in run 
B3-hr is more centrally-peaked than its low-resolution 
counterpart. Since the ACF is simply the Fourier trans- 
form of the power spectrum, this indicates that the high- 
resolution run gives a broader power spectrum with more 
power on the small scales. In addition, comparing vari- 
ous profiles in Figure [3] between runs B3 and B3-hr, we 
see again that the profiles are very similar, with higher 
resolution giving slightly higher turbulent energy. This 



3.2. Energetics and Angular Momentum Transport 

There are two sources of angular momentum transport 
in accretion disks, namely, the radial transport (or re- 
distribution) via turbulence and the vertical transport 
via disk wind. We study them separately in this subsec- 
tion and then briefly discuss about the energetics in our 
simulations. 

3.2.1. Radial Transport of Angular Momentum 

The radial transport of angular momentum is char- 
acterized by the a parameters defined in equation ([7]). 
In Figure we show the vertical profiles of Maxwell 
and Reynolds stresses for all our simulations. Both 
the Maxwell and Reynolds stresses increases monoton- 
ically with net vertical magnetic flux, and approaches 
a few times 0.1po c s f° r Pa — 10 2 - In the mean time, 
the profile gradually changes from centrally peaked for 
large f3 , to a more or less flat profile at f3 = 10 2 . 
In particular, the slope of the wing beyond ±2H be- 
comes more level as /?o decreases. We can further com- 
pare our run B4 with zero net vertical flux stratified 
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Fig. 6. — Vertical profiles of Maxwell stress (solid) and Reynolds 
stress (dashed) for runs B2, B3-hr and B4 (from top to bottom 
panels). In each panel, we show the total stress (bold) as well as 
contribution from turbulent fluctuations (thin). The stresses are 
normalized by po c s- 

MRI simulations ( t atong et al.ll996 HMiller fc Stonell200fl: 
iGuan fc Gammid 120 lit bimon et al.l 120121) . where the 
slope of the wing is further steeper. These results sug- 
gest a general trend that the surface region plays a more 
and more important role in the radial angular momen- 
tum transport as the net vertical magnetic field increases. 
Also, we see that the Maxwell stress is always larger than 
the Reynolds stress, but the Reynolds stress becomes 
more and more important toward disk surface. 

To further decode the mechanism for the radial trans- 
port of angular momentum, we separate the contribu- 
tions from the mean field/mean flow (pv x v' y or Reynolds 

stress and —B x B y for Maxwell stress) and the turbulent 
stress (the rest). As net vertical field increases, the tur- 
bulent component increases steadily with net field and 
saturates for j3 < 10 3 in a way similar to the turbu- 
lent energy density profile in Figure [3l There are two 
intriguing features Figure [5] First, in the disk surface 
region (near our vertical boundary), the stress is dom- 
inated by the mean field/mean flow components rather 



than turbulence, and the contribution from the mean 
field/mean flow components rapidly increases as net field 
increases. In run B4, the mean field component domi- 
nates the wings of the Maxwell stress, while in runs B3-hr 
(and B3), the mean field component dominates the entire 
disk. Second, for sufficiently strong background magnetic 
field (3o < 10 2 , the turbulent stress become completely 
unimportant at all locations, and the turbulent Maxwell 
stress can even become negative. Interestingly, the cen- 
tral region where the turbulent Maxwell stress is positive 
coincides with the bump in the density fluctuation (Fig- 
ure g]). 

The decline of the turbulent component and rise of 
the mean field/mean flow component od the stress as 
one increases net vertical flux strongly contrast with 
conventional unstratified shearing-box simulations. In 
those simulations, the initial net toroidal magnetic field 
is mostly set to zero. Although the net toroidal mag- 
netic flux is a not conserved quantity, it generally stays 
very close to zero for the duration of the simulations. 
With stratification, the large-scale toroidal field evolves 
substantially and for /3q < 10 3 , it well exceeds equiparti- 
tion strength (see Figures [71 and [T2")) . Although the linear 
growth of the MRI is not directly coupled to the toroidal 
field, the non-linear saturation certainly depends on it 
(jHawlev et al.l [1995T) . Therefore, it is not appropriate 
to compare unstratified MRI simulations with stratified 
simulations with the same (3q unless the unstratified sim- 
ulation contains a substantial net toroidal magnetic field. 

Integrating the stresses one obtain the a parameter de- 
fined in (J7J) , and the results are listed in Table [5J The 
a parameter increases monotonically with net magnetic 
flux, from a total of about (a) w 0.075 for (3o = 10 4 to 
(a) w 1.0 for Pq = 100, The ratio of (ctMax) to (aR oy ) in- 
creases from about 4 for j3o = 10 4 , which is si milar to that 
in ze ro net vertical flux simulations (e.g., iDavis et al.l 
l2010f ). to about 7 for /3 = 10 2 , which is mainly due to 
an increasing contribution from the mean (rather than 
turbulent) magnetic field. 

We note that in run B2, the Maxwell stress only de- 
creases very slowly with height, while the Reynolds stress 
even increases with height, making the volume integrated 
stress limited by the vertical extent of our simulation box. 
However, we argue that increasing the box size does not 
necessarily improve the situation. In this case, the radial 
angular momentum transport is dominated by the mean 
field and mean flow. As we shall discuss in Section0] the 
strength of the large-scale mean field at disk surface is 
intimately connected to the global condition of the disk 
and can not be determined in the local shearing-box ap- 
proch. Therefore, uncertainty still remains as one uses a 
taller box. The measured (a) value for run B2 should be 
only taken as a reference. 

3.2.2. Vertical Transport of Angular Momentum 

Angular momentum can also be extracted from the 
disk directly in the vertical direction through outflow and 
magnetic fields, which is determined by the zcf) compo- 
nents of the Reynolds and Maxwell stresses exerted at 
vertical boundaries 

= (pVyV z )\ hot .top , 
rpMlLX f RDM 

J-zcf, — \ — -Dy-t^JIbot.top , 
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TABLE 2 

Energetics and Angular Momentum Transport. 



Run 


("Max) 


( Q Rcy) 








E K 


Ep 


rh w 


B2 


0.92 


0.12 


0.061 


0.048 


3.45 


0.89 


1.26 


0.0944 


B3 


0.37 


0.077 


0.0085 


0.0041 


1.40 


0.079 


0.15 


0.0138 


B3-hr 


0.47 


0.086 


0.010 


0.0046 


1.76 


0.086 


0.21 


0.0160 


B4 


0.061 


0.014 


0.00072 


0.00036 


0.23 


0.0065 


0.010 


0.0013 



Note: All in natural units, with po = c s = Q = 1. 

where in practice subscripts bot and top mean that we 
measure these quantities at the last layer of the top and 
bottom of the computational domain. The torque ex- 
erted by the wind stresses can be obtained by simply mul- 
tiplying (T^ y ) and (T^ ax ) by the radius R. For steady 
state accretion, and include the contributions from radial 
and vertical transports, the accretion rate can be written 
as 



M = AL 



AnR r 



'-T z , 



(11) 



where M r and M z represent contributions to the accre- 
tion rate from radial and vertical transport respectively. 
Qualitatively, one may replace the vertical integral J dz 
by the disk scale height H . We see that given the same 
level of stress T r ^ and T 2( ^, vertical transport is more 
efficient than radial transport by a factor of about R/H. 

Being a local shearing-box simulation, however, R is 
unspecified, and moreover, the location of the central 
object is unspecified (can either at inner or outer x di- 
rection), making the sign of R ambiguous. A physical 
disk wind requires that the sign of T Z( p on the two ver- 
tical boundaries be the opposite and are steady, and is 
closely related to the symmetry of the wind solution as 
we will discuss in detail in Section r4.3l In this subsection, 
we only consider the absolute value of T Z( p and average 
it over the two vertical boundaries (i.e., assuming the 
flow structure and magnetic configuration has the phys- 
ical geometry) 5 , hence the relative importance between 
radial and vertical transport can be rewritten as 



M r 
M~ z 



2ttH 



4 R\T%°\/Poc 



(12) 



We list the values of T^ oy and T™ ax from each sim- 
ulation run in Table [5J with both normalized by po c s- 
The wind stresses appear to depend more sensitively on 
Po, and increase by about two orders of magnitude from 
[3 = 10 4 to P — 100. This may suggest the importance 
of wind transport relatively to turbulent transport in- 
creases with net vertical field. Taking H/R = 0.1, we 
find that wind transport, if present, should be small (less 
than 25%) compared with radial transpo rt for /3p = 10 4 , 
which is consistent with the results by iFromang et al.l 
(|2012f) . It would become more or less comparable to ra- 
dial transport for /3q = 10 3 , while it would dominate 

5 We will discuss in Scction l4.3l that the outflow is unlikely to be 
directly connected to an ordered disk wind mainly due to geometric 
reasons. Let us set aside the issue with outflow geometry here for 
the purpose of discussion. 



radial transport for (3q = 10 2 . In addition, the contri- 
butions from Maxwell T]^ ax and Reynolds T^ ey compo- 
nents are roughly 2 : 1 in runs B3 and B4, while the 
Reynolds component becomes more important for run 
B2. However again, as we shall discuss in Section @J the 
wind stress can only determined from global approach, 
while our measured (T Z( f,) values should also be taken as 
a reference. 

3.2.3. Energetics 

Using an isothermal equation of state, total energy is 
not conserved in the simulation. However, it is impor- 
tant to examine the energy budget for consistency such 
that the work done by the Keplcrian shear (i.e., from 
radial shearing-box boundaries) exceeds the energy loss 
from the vertical boundaries, with the rest of the energy 
presumably escape from the disk in the form of radiation. 

The work done by the shearing-box radial boundaries 
per unit area per unit time is given by 



Psh = - 



1 f 3 3 
— — / -(B x B y }flL x dydz = -ECjO(a Ma 

->X-L>y J Z I 



(13) 

where S = J {p{z))dz is the column density. 
The energy loss from vertical boundaries is given by 

E = E K +E P = pf^-+cl\v-(vxB)xB n 

- \ ^ / J top+bot 

(14) 

where n is the unit vector pointing away from the disk 
in the vertical direction, and we sum over contributions 
from the top and bottom vertical boundaries. The first 
terms in the bracket represent the kinetic energy loss 
and the PdV work done by the mass outflow, and the 
second term represents energy loss from the Poynting 
flux. We have ignored the energy loss term associated 
with internal energy due to the mass outflow, since we 
keep feeding mass to the system which balances this term 
exactly This term is comparable to the PdV work term 
which, as we shall see, is always small compared with the 
total E (hence our mass addition do not substantially 
alter the energy balance). 

In Table [2] we also list the values of P s h, Ek and Ep, 
normalized in natural units. We see that P sri is always 
larger than the sum of Ek and Ep, hence energy con- 
servation is not violated. Also, we sec that Ep is al- 
ways larger than Ek by a factor of about 1.5 to 2, hence 
Poynting flux dominates energy loss. In addition, under 
the assumption of isothermal equation of state adopted 
in our simulations, we find that as net vertical magnetic 
field increases, Ek + Ep roughly increases as l//3n, which 
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is much faster than the rate P s h increases. This may sug- 
gest radiative efficiency decreases as a larger fraction of 
work done by the shear is converted into the wind energy 
loss. 

3.3. Dynamo Activities 

It is well known in stratified shearing-box simula- 
tions with zero net vertical magnetic flux that the mean 
toroidal field component exper iences periodic flips on the 
time scale of 1 orbits (e.g. iBrandenburg et ah 119951 : 
iShi et al.l [20101: iDavis et alJ 12011 " ISimon et al. 1 1201217 . 
Since a zero background vertical field does not give any 
instability on its own in the linear theory, the MRI activ- 
ity in these simulations together with the periodic flips 
strongly suggest that a mean-field dynamo is in opera- 
tion. It is generally believed that shear, turbulence and 
large scale azimuthal fi eld are responsible for sustain- 
ing the dynamo cycles (IVishniac fc Brandenburg IT997t 
lYousef et al.ll2008t iLesur fe Qgilvidl2008l: lGresselll2oToh . 
yet the detailed mechanism remains unsettled. 

In the presence of net vertical magnetic field, the MRI 
can be sustained from the linear modes hence a dynamo 
is not necessary to explain the turbulence. With very 
weak vertical n et magnetic flux, dynamo behavior was 
also o bser ved in [Suzuki fe Inutsukal (|2009f ) , ISuzuki et al.l 
([20101) and lFromang et alJ^OlC whose simulations cor- 
respond to /3q = 10 4 — 10 6 , although it was not studied in 
detail. Our simulations allow us to systematically study 
the evolution of the MRI dynamo with net vertical mag- 
netic flux. In Figure[7j we show the space-time plot of the 
horizontally averaged azimuthal magnetic field for all our 
simulation runs. It is obvious that for (5q = 10 4 , the mean 
azimuthal field still undergoes cyclic oscillations between 
positive and negative signs, while for /3q = 100, the dy- 
namo behavior disappears, where the mean azimuthal 
field in the entire simulation box is completely domi- 
nated by one sign at all times. The case with (3o = 10 3 is 
somewhat marginal, where the dynamo-like behavior is 
present for part of the time while for the rest of the time 
the mean azimuthal field is predominantly one sign. 

In zero net vertical flux stratified shearing -box simu- 
lations (see Figure 15 of ISimon et al.l [20121 for a most 
clear rendering), the dynamo pattern is highly repeat- 
able in the space-time plot, also known as the "butterfly 
diagram" , with a well-defined of period about 10 orbits. 
The dynamo cycles in our run B4 with (3q = 10 4 , how- 
e ver, is highly irregular . Such behavior is also observed 
in lFromang et a l. (2012). Further including the marginal 
case of run B3 with (3 = 10 3 , we see that there is a 
systematic trend that as net vertical field increases, the 
dynamo cycle weakens by becoming more irregular with 
less periodicity. The flipping of mean azimuthal field be- 
comes more sporadic, and the mean time interval for field 
flipping also becomes longer. As (3q reaches below 10 3 , 
the flipping time interval virtually becomes infinity, and 
the dynamo is completely quenched. 

We note that the most unstable linear MRI mode is 
properly resolved in run B2, which does not show dy- 
namo behavior, but may not be well resolved in runs B3 
and B4, which show dynamo activities. Hence question 
arises on whether the dynamo activity depends on the 
proper resolution of the linear MRI modes. Our run B3- 



hr, which is the same as run B3 but properly resolves the 
most unstable MRI mode, is designed to clarify this po- 
tential ambiguity. We see that the space-time pattern for 
the mean azimuthal field is very similar in the two runs: 
dynamo activity appears only part of the time. This 
comparison further justifies that the dynamo behaviors 
observed in our simulations are real. 

It is natural to ask about what physical effects control 
the dynamo activities. Phcnomcnologically, we see that 
in the presence of the dynamo, the disk is separated into 
two distinct regions, namely the region near the mid- 
plane with relatively weak magnetic field where the dy- 
namo appears to be developed, and a highly magnetized 
region outside. In Figure [71 we have overplotted white 
contours that mark the total plasma f3 = 1 (i.e., total 
magnetic pressure equals the gas pressure). Remarkably, 
these white contours perfectly separate the two regions 
for both runs B3 and B4, while for run B2, the entire 
disk is magnetically dominated hence there is no contour 
at all times. It is intriguing to notice that the dynamo is 
present whenever the magnetic field strength in the disk 
midplane region is below equipartition strength. 

Meanwhile, the criterion of f3 = 1 also applies to 
magnetic buoyancy (i .e., the undulatory Parker in- 
stability, iParkerl |1966|) as discussed i n a n umber of 
zero net-flux simulations dBlaes et al.l [20071 IShi et alJ 
[20101: iGuan fc Oammiei 120111 : ISimon e tall l201~lh . With 
isothermal equation of state, and assuming zero 
mean vertical field, the fluid is buoyantly unstable 
if the magnetic energy density decreases with height 
(jGuan fe Gammiell2011l) . It turns out that the magnetic 
field strength tends to be constant in the gas pressure 
dominated disk midplane regions, presumably due to effi- 
cient turbulent mixing, while it falls off with height in the 
magnetically dominated corona. Correspondingly, for re- 
gions with /3 < 1, the fluid is unstable due to magnetic 
buoyancy, while for regions with j3 > 1 , the fluid is buoy- 
antly (marginally) stable and dynamo activities produce 
cyclic alternations of the mean magnetic field. Although 
the analysis of the Parker instability in previous works 
assumes zero mean vertical field, it continues appear to 
be valid in our simulations with net vertical field, which 
is likely because the net vertical field in our simulations 
is much weaker (by at least a factor of 10) than the mean 
azimuthal field. 

The strong azimuthal magnetic field in our simulations 
with /3o < 10 3 is t o some extent simi l ar to the simula- 
tions performed bv I Johansen fe Levinl ()200 8). They ini- 
tiated their simulations by a pure azimuthal field with 
equipartition strength (with zero net vertical magnetic 
flux), and found that the interplay between the Parker 
instability at disk surface layer and the MRI near the 
midplane effectively produces a magnetic dynamo, and 
the azim uthal flux is we l l conf ined to the disk. Simu- 
lations of iGaburov et all (|2012h for the tidal disruption 
of a molecular cloud approximately realize the situation. 
The key difference in our case is that a strong outflow 
is produced in the presence of net vertical magnetic field 
(see Section^]), hence azimuthal flux continuously escape 
the simulation box and has to be continuously generated 
within the disk. Therefore, although Parker instability 
is likely to be present in the entire disk when (3q < 10 3 , 
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Fig. 7. — The time evolution of the horizontally averaged profiles of toroidal magnetic field in all simulation runs, from top to bottom: 
B2, B3, B3-hr and B4. Here we mark t' = as the time when full net vertical flux is added (t = 240ft- 1 for B2, B3 and B4, t = for 
B3-hr). The white contours mark plasma (3 = 1. 



the dynamo mechanism of Johansen & Levin no longer 
operates in the presence of strong net vertical magnetic 
field. On the other hand, magnetic dissipation due to 
the Parker instability is likely to play an important rol e 
on the thermal structure of the disk ([Uzdenskvl I2012D . 
which deserves future exploration. 

4. SIMULATION RESULTS: OUTFLOW 



It has been foun d in ISuzuki fc Inutsukal (|2009l ) and 
ISuzuki et al.l (|2010| ) that the inclusion of a net vertical 
magnetic field leads to strong mass outflow in shearing- 
box MRI simulations, and the rate of the mass outflow is 
approximately proportional to l//3o. The mass outflow 
was interpreted as the launching of a wind from the disk 
and w ould potentially connect to a iBlandfor d fc Pavnd 
(|1982l ) type global disk wind. The connection between 
MRI-driven outflow and global disk wind has been fur- 
ther investigated by iFromang et all ()2012l ). who focused 
on the case with /3 = 10 4 . They argue that the MRI 
is able to launch the magnetocentrifugal wind that is 
strongly time-dependent. S imilar conclusion has been 
made by iLesur et al.l (|2012 f) for the case with (3q ~ 10, 
where there is only one single MRI mode fitted into the 
disk with a physical symmetry for the wind, and sec- 
ondary instabilities are likely to make it time-dependent. 
Very strong outflows are also observed in our simulations. 
However, we argue that the outflow seen in shearing-box 
MRI simulations is unlikely to be directly connected to 
a global disk wind, as we elaborate below. 

4.1. Structure of the Outflow 
4.1.1. Critical Points 

The global disk wind model by IBlandford fe Payne] 
()1982t ) states that when the magnetic field at the surface 



of a thin disk is inclined relative to the disk normal for 
more than 30°, outflowing gas can be accelerated cen- 
trifugally along field lines, extracting angular momen- 
tum from the disk, and is gradually collimated by the 
magnetic hoop stress. A crucial ingredient of this pic- 
ture is the wind launching/mass loading from the disk, 
which is intrinsically connected to the gas dynamics in 
the disks. The wind launching process has been studied 
analytically under the local shearing-box approximation 
(iWardle fc Koenigll Il993t iQgilvie fc Livkll2nol lOgilvid 
l2012h , where all the radial gradients except the shear are 
omitted to reduce the problem to fully one-dimensional 
(ID). The ID solutions are then matched to global wind 
solutions, where the mass loading rate is determined as 
an eigen-value problem. 

A physical wind solution should pass three critical 
points, namely, the slow and fast magnetosonic points, 
and the Alfven point. In the ID (laminar) model, they 
are given by 



1 



1 



n = ^{ci+v\)±^{ci + v\) 



^s u Az ' 

for fast and slow magnetosonic points, and 



±VAz 



(15) 



(16) 



for the Alfven point, where va z = wBf/p is the ver- 
tical component of the Alfven velocity, and va is the 
full Alfven velocity. The requirement that the flow 
passes smoothly through these critical points poses three 
eigen-value problems that would determine three physi- 
cal quantities, such as the mass loading rate and the field 
inclination at disk surface. In the local shearing-box ap- 
proach, the fast magnetosonic point, and sometimes the 
Alfven point, are sufficiently high above the disk where 
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Fig. 8. — Horizontally-averaged vertical velocity profiles for fidu- 
cial resolution simulations B2, B3 and B4 with /3q = 10 2 (solid), 
10 3 (dashed) 10 4 (dash-dotted) respectively. Also marked are 
the locations of Alfven points (solid dots) and slow magnetosonic 
points (circles). The fast magnetosonic points are beyond the sim- 
ulation domain. 

the local approximation no longer applies, and they need 
to be captured in the global models. The lack of one 
or two critical points implies that the solution obtained 
within the computational domain has one or two degrees 
of freedom. 

In Figure [51 wc show the time and horizontally aver- 
aged profiles of the gas vertical velocity as a function 
of height. In all cases, the vertical velocity rapidly in- 
creases towards disk surface, and becomes transsonic or 
supersonic as the flow leaves the simulation box. We also 
calculate the location of the critical points based on equa- 
tions (|T3j) and (Unj) , where the Alfven velocity is based on 
the time average of the absolute value of the horizontally 
averaged the magnetic field. We find that for all cases, 
the fast magnetosonic points are beyond our simulation 
box (the fast magnetosonic speed is about 3c s or larger 
at vertical boundaries), while the Alfven (marked in the 
Figure) and slow magnetosonic points are well contained 
in the box. This fact means that the flow structure is 
not fully determined and has one degree of freedom. 

4.1.2. Conservation Laws 

It is well known that for a laminar magnetized 
wind, the gas follows the magnetic field lines, and the 
flow is characterized by a n umber of conserved quan- 
tities along the streamli nes (|Blandford fc Pavnel Il982t 
iPelletier fc Pudritzlll992[ ). These conservation laws pro- 
vide very useful diagnostics that help us understand 
the mechanism for the launching and acceleration of 
the outflow, as well as the angular momentum trans- 
fer processes. The forms of these conserva tion laws in 
the sh earing-box framework are derived in iLesur et al.l 
(|2012f) . In shearing-box, it was shown that poloidal gas 
streamlines do not necessarily follow exactly the poloidal 
field lines, which brings new terms to the conservation 
law. Nevertheless, the difference is generally small. Since 
strong turbulence are present in all our simulations while 
these conservation laws are derived for a laminar flow, 
they only hold approximately and serve for diagnostic 




Fig. 9. — Top: time averaged vertical profile of mass outflow rate 
normalized to the mean vertical magnetic field (or k, see Equation 
1171 for all our simulation runs. Bottom: time averaged mass loss 
rate from the simulation domain as a function of net vertical flux. 
Thick/thin symbols correspond to fiducial/high resolution simula- 
tions (B3-hr/B3). 
purpose. 

The starting point is mass conservation, where 

pv z = const = kB z , (17) 

where overline indicates horizontal average, and k is con- 
stant. We show the profile of ~pvz~/B z (or k) for all our 
simulation runs in the top panel of Figure [9] Because we 
constantly add mass to the simulation box to maintain 
steady state, the vertical profile of pwl does not become 
flat until beyond z ~ ±3iif . There is some asymmetry in 
the high-resolution run B3-hr, which is most likely due 
to its shorter run time and the lack of statistics. The 
asymptotic values of ~pvl are used to calculate the mass 
loss rate m w from the simulations 

m w = (pwl)ltop - (pwT)lbot , (18) 

and we have listed the value of m w in Table [2J Further 
details about the mass outflow rate will be discussed in 
Section [O] 

As long as k is constant, as valid beyond z ± 3H in 
our simulations, specific angular momentum and energy 
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are co nserved along streamlines. Following iLesur et al.l 
<|2012T ) . the specific angular momentum reads 



Run B2 



Run B4 



K 



(19) 



where C 



fix/ 2 is the fluid part of the specihe 



angular momentum. The partition between C and B v /k, 
then describes the angular momentum exchange between 
gas and magnetic field. 

Energy co nservation is given by the Bernoulli invariant, 
which reads ([Lesur et al.ll2012l without correction terms) 



Ebc 



Ek 

r.2 



Et + Ea, + Eb 



— + 4 iog(p) + 4> - 



(20) 



where E^cr represents the specific energy along a stream- 
line, with the four terms denoting kinetic energy, en- 
thalpy, potential energy and the work done by the mag- 
netic torque, respectively, and 



kB v /p 



(21) 



Note that the full velocity u (rather than v) enters Ek, 
and the potential energy for the shearing-box is 



2 



-Viz 2 



(22) 



We extract the mean flow of the gas and integrate to 
obtain the streamline. Because the disk midplane re- 
gions are generally the most turbulent with n varying 
with height, we do not expect conservation laws to hold 
in these regions region and only trace streamlines from 
z = 3H to z = 6H, setting x = at z = 3H (note that 
by construction the conservation laws do not depend on 
the choice of the zero point of x). This treatment cov- 
ers the most interesting surface regions of the disk where 
the launching and acceleration of the outflow take place. 
For illustration purpose, we only consider run B2 and B4. 
For run B2, the large-scale flow and magnetic structures 
are more or less steady (as seen from Figure [7]), hence we 
take the long-term average to obtain the desired physi- 
cal quantities along streamlines. For run B4, due to the 
dynamo activities, we only pick up a short period be- 
tween t' = 560fi _1 and t' = QOOil^ 1 , where the magnetic 
structure in the upper side of the disk is quasi-steady, as 
judged from Figures [71 and [TT1 

In Figure [10] we show the vertical profiles of individual 
terms in the specific energy and angular momentum for 
streamlines obtained from runs B2 and B4. We see from 
the black dashed lines that energy and angular momen- 
tum are approximately conserved beyond about 4H for 
both cases. The rise of the blue lines in the energy plots 
indicate flow acceleration. This is mostly compensated 
by the reduction of potential energy (red) and work done 
by the magnetic torque (green). Similarly, in the bot- 
tom panels of Figure QJjJ we see that the increase of fluid 
angular momentum is compensated by the reduction of 
magnetic torque. Therefore, it becomes clear that the 
acceleration is magnetocentrifugal in nature. Note that 
the role played _by centrifugal potential 4> and the mag- 
netic torque —B y v*/n can be exchanged by shifting the 
zero point of x, hence they represent the same effect. 




Fig. 10. — Energy (top) and angular momentum (bottom) along 
a streamline in the disk surface from z = 3H to 6H for run B2 
(left, with long time average) and B4 (right, averaged between t' = 
560^-! and t' = W0Q- 1 ). For the energy plot, the various terms 
in Equation J 20 D arc represented by: E^ eI (black dashed), Ex 
(blue solid) , Et (magenta solid) , Eb (green solid) and <f> (red so lid) . 
For the angular momentum plot, the various terms in Equation Q19t 
are represented by: / (black dashed), C (blue solid), —B y /n (red 
solid). 

For example, if we were to trace the streamlines from 
the disk midplane for the case of run B2, then we would 
found that the acceleration is almost completely due to 
the centrifugal potential, since for Run B2 the poloidal 
field lines are almost always inclined for more than 45° 
relative to disk normal, as can be seen in Figure 1121 

4.2. Mass Loss Rate 

In the bottom panel of Figure[9l we show the total mass 
outflow rate M w measured from the two vertical bound- 
aries for all our simulations, with the numbers given in 
Table O We see that rh w sc ales roughly linearly with 
1 //3 n, which is consisten t with ISuzuki fc Inutsukal ()2009l ) 



and lSuzuki et al.l (|2010D , while extending the range of po 
from 10 4 in their case down to 10 2 . From the definition 
of the Alfven point, the mass loss rate can be expressed 
as 



= 2VP 



\B, 



(23) 
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where \a denotes value at the Alfven point. A linear 
scaling of rh w with l/flo indicates that p oc B z oc l//?o 
at the Alfven point, which is in line with the fact that 
the location of the Alfven points decreases in height as 
net vertical field increases. Moreover, based on the same 
scaling, v z at the Alfven point is expected to be constant, 
which is roughly the case as seen from Figure [5] 

4.2.1. Determining m w 

In Figure HI we also indicate the s c aling relation of 
m w with [3 from ISuzuki k Inutsukal (|2009f ). We find 
that while our measured m w follows the same trend, 
while our proportional coefficient is about a factor of 5 
smaller than theirs. Besides that our vertical box size 
is slightly larger (12H versus 8\/2 11.3iT), this dif- 
ference is mainly due to the different outflow boundary 
conditions implemented in our simulations from theirs. 
In fact, the rate of mass outflow rn w is not a well de- 
termined quantity in local shearing-box type simulations 
because of the addition al deg ree of freedom. Moreover, 
in both lFromang et all (|2012h and lLesur et all (|2012[) . it 
was found that rh w decreases as the vertical size of the 
simulation box increases, which again suggests that rh w 
can not be reliably obtained from shearing-box simula- 
tions. 

Since the vertical gravity increases linearly with z in 
shearing-box approximation, the gas has to overcome 
stronger and stronger potential well in order to flow out 
as box height increases. In reality, the vertical increases 
only as z/(R 2 + z 2 ) 3 / 2 , where R is the distance to the 
central object, hence shearing -box approximation fails at 
z > R. iFromang et al.1 (|2012l ) considered higher order ex- 
pansion of the real potential which would also require the 
radial potential to be modified to the same order. They 
noticed that this procedure would introduce curvature 
and no longer fit into the shearing-box framework. 

The mass loss rate in real systems should be a well de- 
fined quantity, and we speculate that it should be com- 
parable to the mass loss rate in shearing-box simulations 
performed with box size L z sa 2R. Suppose the mass 
outflow rate scales as 

m w = A Bl 
2p c s (L Z /H)»2P ' 1 ; 

where A is a dimcnsionless constant, the index v de- 
scribes the scaling of rh w on the box height, and we have 
assumed that m w oc l//?o as suggested by the simulations 
discussed in the previous subsection, with Bq denoting 
the background vertical field strength. With a little al- 
gebra, we find 

In particular, if v w 1 (S. Fromang, 2011, private com- 
munication), then our hypothesis (L z = 2R) yields 
rn w sa ABq/2RH. This means that the mass loss rate 
can be found once the coefficient A is known. Moreover, 
m w is solely determined by the background vertical field 
strength and does not depend on the disk thickness or 
the disk surface density. On the other hand, if v deviates 
from 1, then one would expect m w oc (H/R) v ~ l . 



In sum, although rh w is not well determined from 
shearing-box simulations, we argue that by carefully 
studying the dependence of m w on the height of the sim- 
ulation box, it is possible to obtain a physical estimate 
of the mass loss rate. 

4.2.2. Evolutionary Scenarios 

While we have demonstrated the correlation rh w oc 
l//?o, we have not discussed the range that it is applica- 
ble. Reading from Table El we see that rh w for /3q = 100 
is somewhat smaller than the expected linear correla- 
tion, which is suggestive of saturation. Indeed , the rate 
of mass outflow measured in lLesur et al.l (|2012f) for a sim- 
ilar numerical setup (their run 1Dz6) with [3 ~ 10 gives 
m w = 0.234 (multiplied by 2 to account for two sur- 
faces), which is only a factor of 2.5 times higher than 
our case with /3 = 100 {rh w = 0.0944). With fur- 
ther stronger m agnetic field th at is too strong for the 
MRI to operate, lOgilviei (|2012l ) showed that the rate of 
mass outflow should rapidly decrease once /?o drops be- 
low 1. Therefore, the scaling relation m w oc l//3o holds 
up to /3q > 100, beyond which the mass loss rate slows 
down and will eventually fall off rapidly for strong super- 
equipartition field field. 

These results, all together, suggest an evolutionary sce- 
nario. Assuming no magnetic flux evolution, the disk 
would lose mass at constant rate (assuming v = 1) un- 
til the midp l ane n et vertical field exceeds equipartition. 
iLesur et al.l (|2012h explored such a scenario in their ID 
simulations and confirmed the steady decrease of mass 
loss rate with increasing magnetization in the strong field 
regime. A strongly magnetized thin accretion disk may 
gradually evolve into such a stage, where the rate of mass 
loss is self-regulated so that it loses mass faster as more 
mass is fed into the disk (increasing /3q, but still /3q < 1), 
while loses mass slower if the mass loss is not com- 
pensated (decreasing (3 ). T he disk structure would b e 
similar to a jet-emitting disk (jCombet fc Ferreirall2008l ) . 
which is tenuous with super equipartition field. Alter- 
natively, if the evolution starts from a strongly magne- 
tized thick disk, the mass loss rate can be so large that 
the disk will be depleted in a few orbits. For example, 
without feeding mass to our run B2, the mass loss time 
scale is only about 5 orbits. Such rapid mass loss may 
lead to catastrophic disruption of the disk, such as in 
core-collapse supernovae where rapidly rotating progen- 
itor core is threaded by extremely strong magnetic fields 
([Burrows et al.ll2007t ). 

4.3. Fate of the Outflow 

In this subsection we discuss whether the outflow 
launched from shearing-box MRI simulations can be con- 
nected to a global Blandford-Payne type disk wind. The 
key to this question lies in the symmetry. We note that 
due to the neglect of curvature, shearing-box approxi- 
mation does not specify which side of the radial domain 
the central object is located. For a physical disk wind, 
symmetry requires that the outflow from the top and bot- 
tom sides of the disk should be inclined towards the same 
radial and azimuthal directions (which would be consid- 
ered as radially outward and trailing) , and remains so at 
all times. Since gas flows along magnetic field lines, the 
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Fig. 11. — Same as Figure [7] but for run B4 only, and color represents the horizontally averaged radial magnetic field. 

field lines at the top and bottom sides tend to incline 
toward different directions (as can be tracked in Figures 
EandHTJ). 

Symmetry is also directly related to the vertical an- 
gular momentum transport by the outflow, as we dis- 
cussed in Section 13.21 With the wrong symmetry, the 
wind stress (|10[) at the top and bottom sides of the disk 
have the same sign, and cancel each other, hence the net 
vertical angular momentum transport vanishes. 

One might attribute this "wrong-symmetry" problem 
to the one extra degree of freedom in shearing-box sim- 
ulations, which could possibly be avoided by perform- 
ing global simulations. While it is certainly important 
to explore the sa me probl e m in global geometry, it has 
been suggested by lOgilvid (|2012f ) that global effects can 
be effectively taken into account by applying some ad- 
ditional constraints on the vertical boundaries. Ogilvie 
constructed a ID model for launching disk wind / jet in 
the shearing-box approximation, where a much stronger 
net vertical field /3 < 1 was assumed so as to suppress 
the MRI. He assigned the two free parameters at the ver- 
tical boundaries to be the strength of the radial and az- 
imuthal magnetic fields (he has two parameters because 
the Alfven point turns out to be beyond his computa- 
tional domain), and was able to successfully construct 
laminar wind solutions with the desired symmetry for 
the given constraints. 

In view of this possibility, we further modify our ver- 
tical boundary conditions so that the mean radial (or 
azimuthal) magnetic field in the last four grid zones at 
the vertical boundary are set to some fixed value, which 
is done by first calculating the mean radial (or azimuthal) 
field in these zones, and then subtracting a constant ra- 
dial (or azimuthal) field uniformly in all these cells to 
make the mean field equal to the desired value with- 
out violating the divergence-free condition. As we see 
in Figure 1121 the mean radial field at vertical bound- 
aries is roughly the same as the mean vertical field in 
strength. Therefore, we demand that the field incline by 
45° at vertical boundaries in the x—z plane, and require 
that the mean radial fields at the two vertical boundaries 
have opposite signs which conforms to the desired sym- 
metry. However, after performing various experiments on 
smaller domain test runs with parameters similar to runs 
B2, B3 and B4, we find that such a treatment does not 
seem to have a physical effect on the profile of the mean 
flow and mean field, except introducing a strong current 
sheet at the vertical boundaries. Even though such addi- 



FlG. 12. — Time-averaged vertical profiles of magnetic field (left) 
and velocity (right) for run B2 with /So = 100. All in code units. 

mean magnetic field at the top and bottom sides of the 
disk should be bent to the same direction. 

The first difficulty for outflow-wind connection is the 
magnetic dynamo, applicable to the cases with (3q > 10 3 
as we discussed in Secion |3"31 The constant flipping of 
the mean azimuthal (as well as radial, sec Figurc fTTj) field 
lines means that the direction of the outflow would con- 
stantly swap between radially inward and outward, either 
of the directions would be invalid for a global disk wind. 

10 4 are in many as- 



Although our simulations with /3q 

pects v ery similar to simulations bv ISuzuki fc Inutsukal 
(|2009f) . ISuzuki et al.l (|2010l) . and lFromang et al.l (|2012T T 

they argue that such phenomenon leads to the time vari- 
ability of the wind. However, we stress here that the 
simple fact of cyclic sign change of the large-scale field is 
already inconsistent with global wind geometry. 

The second difficulty is that even there were no dy- 
namo, such as the case for (3o = 10 2 , the natural sym- 
metry from the simulation is unphysical for it to be con- 
nected to a global disk wind. In Figure [12] we show the 
vertical profiles of the magnetic field and velocity fields 
for run B2. Clearly, the outflows at the top and bottom 
sides of the disk are inclined toward opposite directions, 
which is inconsistent with a global disk wind. Even in 
the presence of the dynamo, there is a good chance that 
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tional constraints at vertical boundaries can in principle 
affect the gas dynamics around disk midplane via fast 
magnetosonic waves, it appears that causality is directed 
from midplane to surface, as the "butterfly" pattern rises 
from the midplane consequent of buoyancy, and such ad- 
ditional constraints beyond the Alfven point can hardly 
affect the dynamics in the midplane. There can be sev- 
eral reasons for the discre pancy be t ween our test results 
and the ID simulations bv lQgilviel ([2012). In particular, 
the Alfven points are not contained in Ogilvie's simula- 
tions hence boundary conditions have direct influence on 
the structure of the solution. Moreover, the reflection 
boundary condition at midplane in Ogilvie's study has 
already enforced the physical symmetry of the solution. 
It is unclear whether that solution is stable if one relaxes 
this constraint and study a full disk. 

In sum, we conclude that whether or not the MRI tur- 
bulence is accompanied by dynamo activities, the outflow 
from the shearing-box MRI simulations is unlikely to be 
directly connected to a global disk wind. The fate of the 
outflow remains to be explored using global simulations. 

5. DISCUSSION AND CONCLUSIONS 
5.1. Summary 

In this paper, we have successfully performed local 
stratified shearing-box simulations of the MRI that in- 
clude a strong vertical magnetic flux with midplane 
plasma f3 of the net vertical field ranging from /?o = 10 4 
to 10 2 , a regime that has not been explored while is very 
likely relevant to accretion disks in many astrophysical 
systems. Such a magnetic configuration gives rise to very 
vigorous MRI turbulence, and simultaneously launches 
an outflow. We studied the properties of the MRI turbu- 
lence as well as the disk outflow in detail and our major 
findings are summarized below. 

For the properties of the MRI turbulence, we find 

• With relatively weak net vertical field of /?o > 10 3 , 
the disk consists of a gas pressure dominated disk 
midplane and a magnetic dominated disk corona 
(P ~ 0.1 — 1). By contrast, the entire disk be- 
comes magnetically dominated when /3q < 10 3 . 
The strong magnetic support substantially modi- 
fies the vertical structure of the disk and makes it 
substantially thicker. 

• Turbulent magnetic and kinetic energies increase 
monotonically with net vertical magnetic flux, and 
saturate when (3 < 10 3 . The MRI also generates 
large scale toroidal magnetic field whose strength 
increases monotonically with net vertical field, and 
dominates the total magnetic energy over contri- 
bution from turbulence for (3 < 10 3 . 

• The Shakura-Sunyacv parameter a increases 
monotonically with net vertical magnetic flux, 
ranging from a ~ 0.08 at (3 = 10 4 to a > 1.0 
for (3 = 10 2 . Maxwell stress dominates Reynolds 
stress by a factor of 4-7. For weak net vertical flux 
of po > 10 3 , turbulent fluctuations dominates the 
contribution to a, while for (3q < 10 3 , radial trans- 
port of angular momentum is dominated by the 
large scale fields. 



• The MRI dynamo that generates cyclic flips of the 
mean toroidal field persists in the presence of weak 
net vertical magnetic flux, but becomes more spo- 
radic with less periodicity as net flux increases. 
The dynamo is completely suppressed when the 
net flux is strong with /3q < 10 3 , where the mean 
toroidal field exceeds cquipartition strength and 
never flips. 

Our results demonstrate the crucial dependence of 
the behavior of the MRI turbulence on the net ver- 
tical magnetic flux. In particular, there is a critical 
net vertical magnetic flux of /3q « 10 3 near which 
many aspects of the MRI turbulence change qualita- 
tively. Additionally, more careful convergence study is 
still needed, and deeper understanding about the prop- 
erties with MRI turbulence would furth er benefit from 
the s t udy of Prantl number dependence (iFromang et all 
I2007t lLongaretti fc Lesurll2010t iLesur et al.ll2012f h 

For the properties of the disk outflow, we find 

• The slow magnetosonic point and the Alfven point 
are always contained in our simulation box. The lo- 
cation of these critical points shifts towards the disk 
midplane as net vertical magnetic flux increases. 
The launching and acceleration of the outflow are 
due to the magnetocentrifugal mechanism, where 
surface magnetic field is sufficiently inclined. 

• There is a robust trend that the outflow mass loss 
rate rh w increases with increasing net vertical field, 
and tends saturate at 0q < 10 2 . Although m w can 
not be well determined from shearing-box simula- 
tions, its exact value in real systems is likely to be 
much smaller than the measured from local simu- 
lations and may depend on the aspect ratio H/R. 

• The outflow from the MRI simulations is unlikely 
to be directly connected to a global disk wind for 
geometric reasons. For j3o > 10 3 , the large scale ra- 
dial and azimuthal field lines are constantly flipped 
due to the dynamo activities without a permanent 
bending direction. For /?o < 10 3 , the large scale 
magnetic fields do not change sign across the mid- 
plane, hence the field lines at the two sides of the 
disk bend to opposite directions, inconsistent with 
a global wind geometry. 

• Angular momentum transport by the disk outflow 
is likely to play a minor role compared with that 
by the MRI turbulence for relatively weak net ver- 
tical magnetic flux. With stronger net vertical flux 
(A) ^ 10 ), the symmetry of the flow (as the dy- 
namo is suppressed) makes net vertical angular mo- 
mentum transport be zero, but gives oppositely di- 
rected inflow/outflow with substantial mass flux. 

It has been reported that rh w from shearing-box MRI 
simulations decreases with increasing height of the simu- 
lation domain, which is related to the intrinsic limitations 
of the shearing-box framework. We propose that a care- 
ful study of the height dependence can help determining 
the true value of rh w in real systems. Our results leave an 
open question on the fate of the outflow, which can only 
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be reliably explored using global simulations and will be 
discussed further in Section [5.31 

5.2. Implications for Global Disk Evolution 

We note that shearing-box simulations with zero net 
vertical flux tend to gi ve a ~ 0.01 — 0.02 with confirmed 
numerical convergence (|Shi et alJ20T?llDavis et alEoifl 
ISimon et al.l 120121) . The vast dominance of simulations 
of this type has, to some extent, implicitly created a 
misconception that the a value resulting from the MRI 
turbulence if of the order 0.01. On the observational 
side, the value of a can be estimated from the viscous 
timescale of transient systems such as dwarf novae, X- 
ray transient and FU Ori burs ts, where it was found to 
be of the order 0.1 o r larger (|Smakl [l999t iDubus et al.l 
20011: IZhu et al.ll2007f) . Such apparent discrepancy leads 
King et all (|2007l ) to question about the efficiency of the 
MRI in driving disk accretion and evolution. Our simula- 
tions demonstrate that the discrepancy can be naturally 
resolved if the disks are threaded by some large scale 
poloidal magnetic fields. 

The strong dependence of a on /3q immediately in- 
dicates that the evolution of accretion disks strongly 
depends on the distribution of poloidal magnetic flux 
through the disk. On the other hand, because of flux 
freezing in ideal MHD, the distribution of magnetic flux 
is intimately connected to disk evolution (as well as 
turbulent diffusion and rcconncction due to the MRI). 
Therefore, the mutual dependence of disk evolution and 
magnetic flux evolution makes the problem of disk ac- 
cretion intrinsically non-local, hence a full understanding 
of accretion disk structure and evolution requires global 
simulations. 

Most global simulations to date adopt a magnetic field 
geometry that is analogous to zero net-flux shearing- 
box: the simulations are initialized either with sin- 
gle/m ultiple poloidal magnetic loops for a thick disk 
torus (|Hawleyll2000T : iDe Villiers et al.ll2003HPenna et al.l 
120101) . or with pure toroidal magnetic field in a power- 
law disk (iFromang k NelsoiJl2006t iBeckwith et alj|2011t 
iFlock et al.ll2011l) . These simulations generally show val- 
ues of a that are slightly higher than 0.01 (~ 0.025). In 
these simulations, the MRI generates net vertical mag- 
netic flux in the local patches of the disk connected by 
coronal loops and the local stress correlates with the ne t 
vertical flux (|Tout k Pringldll992tlSorathia et al.ll2010D . 
Combined with our study, the local net vertical flux in 
these global simulations is sufficient to account for the 
slightly higher total stress a than that in zero net flux 
shearing-box simulations. However, the distribution of 
local net vertical magnetic flux falls off quickly with 
stronger net flux, hence the total value of a can not be 
enhanced substantially. 

Initial magnetic field geometry plays an important role 
in global simulations of the MRI, although the situa- 
tion with large scale poloidal fields threading through 
the disk have rarely been explored. Our simulations pro- 
vide strong motivation for exploring the consequence of 
large scale poloidal field in global simulations. Such sim- 
ulations will not only expand our knowledge the MRI, 
but will also address the fate of the outflow, which we 
discuss in the next section. 



5.3. Disk Wind Launching in Global Simulations? 

It has been recognized that strong net vertical field 
that approaches equipartition strength at the disk 
midpl ane is necessary for l aunching a steady-state 
wind (IWardle k Koenigj||1993HFerreira k Pelletiedll995t 
iQgilvie k Livioll2001l) . T he requirement may be relaxed 
to lower magnetizations (jMurphv et al.ll2~010l ). while for 
magnetization that is too l ow the outflow is eithe r sup- 
pressed or highly unsteady ()Tzeferacos et al.ll2009D . Nu- 
merous global simulations have studied launc hing of disk 
wind in the context of pr otost ellar disks (e.g..lKato et al.l 
l200llCasse k Keppensl 12001 iZanni et al.1 120071) 1. How- 
ever, as noted in Scction[T] the use of artificial diffusion in 
these simulations do not properly reflect the disk micro- 
physics. For strong net vertical field that would suppress 
the MRI, the artificial diffusion has unknown origin (and 
the non-ideal MHD effects can not be properly treated 
as artificial diffusion. iBai k S tone 2012|). For weaker net 
vertical field, the MRI is the presumed source of artifi- 
cial diffusion. However, our evidence against direct wind 
launching all originate from the microphysical effects of 
the MRI (such as the MRI dynamo) that can not be 
captured in these simulations. 

Another family of global simulations focus on non- 
radiative accretion disks around black holes (BH) . These 
simulations cither use psudo-Newtonian gravitational po- 
tential for the BH or adopt a general-relativistic pre- 
scription. Despite the non-Kcplcrian nature in the vicin- 
ity of the BHs, the physics of the MRI that is not too 
close to the innermost stable orbit is similar. Most of 
these simulations adopt similar poloidal loop field ge- 
ometry and focus on jet launching which deals with the 
innermost part of the accretion flow, a few of them 
do h ave included lar g e scal e poloidal fields. In partic- 
ular, IBeckwith et al.l ((2009) initialized the simulations 
with pure vertical magnetic fields with /3q ~ 10 2 . How- 
ever, upon saturation the launching of a Blandford-Payne 
type wind was not observed, which was attributed to a 
global "coronal mechanism" that governs the magnetic 
flux evolution. Also, their simulations have relatively 
small radial domain and v e ry thick disk. More rece ntly, 
Tchekhovskov etaL[j20 1 If ) . iMcKinney et al.1 (|2012f ) and 
Naravan et al.l ( 20121 ) studied the accretion flow with 
large scale initial poloidal field in the inner disk with 
/3o ~ 10 2 which leads to a "magnetically arrested" or 
"magnetically chocked" accretion flow due to the accu- 
mulation of magnetic flux in the vicinity of the BH. The 
MRI in the main disk was properly resolved in these 
simulations and there is also strong outflow from the 
disk with rate of the outflow increasing with radius, 
which is typical for adv ection dominated accretion flow 
([Stone fc Pringld 120011 ) . Nevertheless, given the com- 
plexity of the flow structure, the origin and nature of the 
outflow is still debated (see also lYuan et al.|[2012f ). 

Together, it appears that Blandford-Payne type disk 
wind has not been unambiguously identified in global 
simulations with properly resolved MRI turbulence, 
which is at least in line with our observations that wind 
launching from MRI-turbulent disks violates geometric 
constraints. More explorations on the magnetic field ge- 
ometry (especially the field geometry with large scale 
poloidal fields) in global simulations and comparisons 
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with local shearing-box simulations are essential to ascer- 
tain the nature of the outflow, and to reveal the under- 
lying physics governing the evolution of accretion disks. 
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